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ABSTRACT:  Numerical  computations ,  using  a  one -dimensional  hydrodynamic 
code  on  an  IBM  7090,  have  been  made  for  the  detonation  of  1-lb  spheres 
of  pentolite  and  TNT  at  sea  level.  Positive  durations,  overpressure  at 
the  main  shock,  and  space-time  paths  of  the  main  shock,  second  shock, 
and  contact  surface  are  compared  with  experimental  data.  It  is  found 
that  the  calculated  peak  pressures  at  the  main  shock  are  10  per  cent 
too  low,  the  experimental  second  shock  path  lies  about  40  per  cent  above 
the  calculated  one,  and  the  experimental  contact  surface  path  is  about  a 
factor  of  two  greater  in  radius  than  the  calculated  one.  Possible  explana¬ 
tions  for  these  deviations  are  discussed. 
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INTRODUCTION 

One  of  the  most  important  methods  available  for  the  theoretical 
prediction  of  blast  effects  in  air  is  the  solution  of  the  hydrodynamic 
equations  by  finite  difference  methods  using  an  electronic  computer. 

The  present  investigation  is  concerned  with  the  blast  wave  and  other 
phenomena  arising  from  the  detonation  of  small  spherical  charges  at  sea 
level  and  their  computation  by  numerical  methods.  In  particular,  we 
present  results  for  the  detonation  of  1-lb  spheres  of  pentolite  and  TNT 
at  sea  level  and  we  compare  these  results  with  available  experimental 
data. 
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THE  NUMERICAL  COMPUTATIONS 

The  hydrodynamic  code  used  in  these  computations  utilizes  the 
Von  Neumann-Richtmyer  artificial -viscosity  technique,1  and  is  fully 
described  elsewhere.8  The  basic  idea  is  to  solve  a  set  of  finite  dif¬ 
ference  equations  (which  approximate  the  partial  differential  equations 
of  gas  dynamics)  by  a  step-by-step  procedure  advancing  in  time.  The 
physical  situation  is  assumed  to  be  spherically  symmetric,  so  that  all 
quantities  are  functions  only  of  the  radius  and  the  time;  the  origin  of 
coordinates  is  chosen  to  be  at  the  center  of  the  explosive.  Hydrodynamic 
variables  are  assigned  values  in  a  set  of  spatial  zones  covering  the  range 
of  interest  at  time  t  =  0,  and  the  flow  evolves  from  these  initial  condi¬ 
tions.  The  origin  of  time  is  taken  at  the  instant  the  detonation  wave 
reaches  the  surface  of  the  explosive,  so  that  the  solid  explosive  has 
been  completely  converted  to  gaseous  detonation  products. 

The  equation  of  state  used  for  the  detonation  products  is  the  ISZK 
equation,3  and  the  initial  conditions  in  the  detonation  products  are 
represented  by  a  spherical  Taylor  wave.3’4  The  equation  of  state  for 
the  air  is  E  3  P/p(Y-l),  where  V  is  a  fit  by  straight  line  segments3  to 
the  tables  of  Hilsenrath  and  Klein,5  and  is  taken  to  be  a  function  of 
internal  energy,  E,  and  density,  p.  The  parameter  v  is  not  equal  to  the 
ratio  of  specific  heats,  except  in  the  low-energy  domain. 

Since  the  artificial-viscosity  mechanism  has  the  tendency  to  smooth 
out  shock  fronts  so  that  they  extend  over  a  small  number  of  spatial  zones , 
it  is  necessary  to  define  how  values  are  to  be  read  at  a  shock.  The 
position  of  the  shock  is  located  in  the  numerical  results  as  the  point 
at  which  the  artificial  viscosity  is  a  maximum,  and  the  pressure  is 
taken  to  be  the  pressure  at  the  first  maximum  (in  pressure)  immediately 
behind  this  point. 
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THE  EQUATION  OF  STATE  FOR  DETONATION  PRODUCTS 

We  have  chosen  to  describe  the  detonation  products  by  the  ISZK 
(Landau,  Stanyukovlch,  Zeldovich,  and  Korapaneets)  equation  of  state, 
which  treats  the  products  by  drawing  an  analogy  between  the  state  of 
the  detonation  products  of  a  condensed  explosive  and  the  crystal  lattice 
of  the  solid  state.  Details  concerning  the  derivation  and  significance 
of  the  LSZK  equation  of  state  are  given  elsewhere,*  and  we  present  here 
only  a  brief  sunmary  of  the  relevant  equations  and  values  for  the  various 
constants.  The  (E,  P,  p)  form  of  the  LSZK  equation  may  be  written: 


P 


(1) 


where  P  =  pressure,  p  *  density,  E  =  internal  energy/unit  mass,  cr  and  V 
are  dimensionless  constants,  and  B  is  a  constant  with  dimensions  (ML-*)1*^ 
(l?T"2).  The  three  undetermined  constants,  a,  y,  and  B,  must  be  evaluated 
by  using  experimental  results.  It  is  sufficient  to  know  the  detonation 
velocity  versus  loading  density  relationship  and  the  value  of  the  heat 
of  detonation  for  the  solid  explosive,  if,  in  addition,  we  postulate 
that  the  detonation  gases  approach  ideal  gas  behavior  at  small  densities. 
The  procedure  for  evaluating  the  constants  has  been  described  elsewhere* 
for  TNT,  and  the  results  are: 


ot  =  2.9412 

|  -  0.53562  (2) 

Y  *  2.78, 


where  Q  =  1018  calories/gram  ■  heat  of  detonation,  and  the  ideal  gas 
value  for  the  ratio  of  specific  heats  of  the  products  at  low  densities 
was  taken  to  be  Yi  =  1.3^*  Using  these  values,  the  spherical  Taylor 
Wave  for  TNT  at  a  loading  density  of  q,  =  I.625  gm/cc  was  calculated 
end  used  for  the  initial  condition  for  the  products  at  t  =  0.  Equation 
(l) ,  with  the  values  of  the  constants  given  in  (2),  is  then  used  for  the 
subsequent  calculation  of  the  state  of  the  detonation  products  for  t  >  0. 

Precisely  the  same  procedure  is  applicable  to  pentolite,  and  has 
been  carried  out  to  yield  the  following  results : 
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a  =  2.9412 


B 


l-V 


|  *  °-”927  (fSr) 

Y  -  2.64, 

where  Q  =  1260  cal/gm  was  used  for  the  heat  of  detonation,  and  the  experi¬ 
mental  curve  for  detonation  velocity  versus  loading  density  was  taken 
from  Cook.9  The  value  of  y*  used  was  1.34  as  for  TNT.  The  initial 
density  of  the  solid  explosive  was  taken  to  be  ^  =  I.67  gm/cc. 


4 


NOLTR  65-57 


COMPARISON  WITH  EXPERIMENTAL  DATA  FOR  PENTOLITE 

A  comprehensive  compilation  and  analysis  of  free-air  blast  data 
from  pentolite  detonations  has  been  produced  by  Goodman,7  the  results 
of  which  will  be  used  in  this  report  to  provide  comparisons  with  the 
theoretical  pentolite  calculations.  The  quantities  we  shall  consider  for 
pentolite  are  the  overpressure-distance  curve  for  the  main  shock,  the 
radius -time  path  of  the  main  shock,  and  the  positive  duration  as  a  func¬ 
tion  of  distance.  Goodman  has  reported  these  quantities  as  single  curves, 
which  are  obtained  by  averaging  and  smoothing  the  results  of  many  dif¬ 
ferent  experimental  determinations.  Since  neither  ambient  pressures  nor 
temperatures  were  usually  reported,  all  measurements  were  arbitrarily 
assumed  to  have  been  made  at  14.7  psi  and  300°K  with  ambient  sound  speed 
1139.4  ft/sec.  The  charge  density  was  assumed  to  be  1.65  gm/cc,  which 
is  close  to  the  value  used  in  the  theoretical  calculation  of  ^  *  I.67 
gm/ cc. 

Figure  1  is  a  plot  of  shock  overpressure  ( total  pressure  minus 
ambient  pressure)  versus  radial  distance.  The  dashed  curve  represents  the 
smoothed  numerical  fit  of  Goodman  to  the  experimental  data,  while  the 
solid  curve  represents  the  calculated  result.  In  the  range  0.2  -  1.5  feet, 
the  behavior  of  the  calculated  curve  is  characterized  by  the  phenomenon 
of  shock  formation,  consisting  of  a  sharp  rise  in  pressure  at  early  times, 
followed  by  an  overshoot  region  in  which  the  calculated  curve  lies  10-20 
per  cent  above  the  experimental  curve.  The  region  is  probably  dominated 
by  purely  computational  effects  which  have  little  to  do  with  the  actual 
physical  phenomena.  Thus,  in  a  real  gas,  the  shock  formation  is  governed 
by  viscous  and  heat  conduction  effects;  however,  in  numerical  calculations 
of  this  type  the  analogous  quantity  is  an  artificial  viscosity  which  is 
quadratic  in  the  space  gradient  of  the  particle  velocity,  rather  than 
linear,  as  is  a  real  viscosity.  The  artificial  viscosity  is  constructed 
to  yield  the  correct  entropy  across  a  shock  but  does  not  reproduce  the 
proper  shape  or  extent  of  the  shock  transition  itself,  nor  does  it 
reproduce  properly  the  physical  phenomenon  of  shock  formation.  Further¬ 
more,  some  difficulty  arises  at  the  interface  between  the  gas  products 
and  the  air,  due  to  the  fact  that  the  zones  on  either  side  of  the  inter¬ 
face  have  widely  different  masses.8  It  would  therefore  seem  unwise  to 
ascribe  too  much  physical  significance  to  the  early  stages  of  the 
calculation. 

At  approximately  1.5  feet  the  calculated  curve  crosses  the  experi¬ 
mental  one;  and  from  then  on  the  two  curves  remain  parallel,  the  calculated 
values  lying  about  10  per  cent  below  the  experimental  values.  The  reason 
for  this  discrepancy  is  not  yet  clear.  Various  possible  explanations 
will  be  discussed  in  a  subsequent  section  of  this  report. 

Figure  2  is  a  radius -time  plot  of  various  aspects  of  interest  in  the 
flow.  The  upper  two  curves  depict  the  arrival  times  for  the  main  shock; 
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the  solid  curve  represents  the  calculated  arrival  times ,  and  the  dashed 
curve  represents  a  smoothed,  numerical  fit  to  the  experimental  data.7 
Also  shown  is  the  calculated  radius -time  plot  for  the  contact  surface 
separating  the  detonation  products  from  the  shocked  air.  Each  of  the 
successive  relative  minima  in  this  curve  represents  a  point  at  which  an 
internal  shock  in  the  detonation  products  passes  into  the  air,  at  the 
same  time  reflecting  another  shock  inward  toward  the  center.  The  first 
of  these  subsidiary  shocks  will  be  considered  in  detail  for  TNT,  in  the 
next  section,  as  experimental  data  exist  for  the  location  of  this  shock 
for  this  explosive. 

The  remaining  curve  is  useful  in  calculating  positive  durations,  and 
represents  the  path  of  that  point  behind  the  main  shock  at  which  the  pres¬ 
sure  is  equal  to  the  ambient  air  pressure  (overpressure  equals  zero).  A 
slight  dip  is  discernible  in  this  curve  at  about  120  p,secs ,  and  is  due  to 
the  fact  that  at  this  time  the  second  shock"' (the  shock  coming  back  through 
the  detonation  gases)  crosses  the  zero  overpressure  curve.  It  is  clear 
from  this  curve  that  the  point  of  zero  overpressure  remains  within  the 
explosion  products  until  800  ^secs,  at  which  time  it  emerges  into  the 
compressed  air  behind  the  main  shock.  At  this  point  a  change  in  the 
slope  of  the  zero-overpressure  curve  occurs.  This  results  in  a  discon¬ 
tinuity  in  the  slope  of  the  positive  duration  versus  time  curve  depicted 
in  Figure  3- 

The  solid  curve  of  Figure  3  is  the  calculated  one,  and  the  dashed 
curve  is  a  smoothed  numerical  fit  to  Goodman's  compiled  pentolite  data. 

The  cusp  does  not  appear  in  the  experimental  curve,  although  the  general 
trend  of  the  two  curves  appears  to  be  the  same.  The  absence  of  the  cusp 
in  the  experimental  data  might  be  explained  by  the  sparsity  of  experi¬ 
mental  points  in  that  particular  region,* **and/or  by  the  fact  that  the 
dashed  curve  is  a  smoothed  numerical  fit  to  experimental  data  which  dis¬ 
play  a  fair  amount  of  scatter. 

In  any  case,  the  agreement  shown  in  Figure  3  is  considered  to  be 
fairly  good  in  view  of  the  fact  that  the  positive  duration  depends  on 
the  conditions  in  the  whole  field  at  successive  times,  and  is  therefore 
subject  to  many  sources  of  perturbation  and  error. 


*  This  shock  originates  at  the  interface  between  the  detonation  gases  and 
the  air.  It  moves  inward  into  the  detonation  gases,  is  reflected  from 
the  center,  and  then  proceeds  outward,  emerging  from  the  detonation 
products  at  about  1750  microseconds. 

♦♦Recent  pressure  measurements,  on  a  large  hemispherical  charge  constructed 
from  32-lb  blocks  of  TNT,  have  been  made  in  this  region.  The  results 
indicate  that  the  cusp  does  exist  for  this  multiton  explosion.  We  are 
indebted  to  Mr.  John  Keefer  of  the  Ballistics  Research  Laboratories  for 
this  unpublished  information.  In  addition,  Adushkin  and  Korotkov15  have 
experimentally  detected  a  cusp  in  the  positive  duration  curve,  although 
they  ascribe  this  effect  to  late  shock  formation. 


6 


NOLTR  65-57 


COMPARISON  WITH  EXPERIMENTAL  DATA  FOR  TNT 

A  comprehensive  analysis  of  experimental  data  for  TNT,  along  the 
lines  of  the  Goodman  report  on  pentolite,  does  not  exist  at  present, 
so  that  we  content  ourselves  with  data  from  particular  TNT  experiments. 

We  shall  make  use  of  experimental  data  from  the  work  of  Fisher9*10  (first 
and  second  shocks),  Weibull11  (first  shock).  Potter  and  Jarvis13  (first 
and  second  shocks),  and  Rudlin13*14  (first  and  second  shocks,  and  con¬ 
tact  surface*).  It  should  be  noted  that  the  loading  densities  were  not 
the  same  for  all  of  these  experiments.  Fisher  used  q,  *1.51  gm/cc, 

Rudlin  used  various  densities  ranging  from  q,  =1.01  gm/cc  to  ^  a 
I.625  gm/cc,  and  the  densities  used  by  Potter  and  Jarvis  were  not  available. 
The  loading  density  used  for  the  theoretical  calculation  was  =  1.625 
gm/cc.  The  calculated  quantities  for  TNT  which  we  shall  compare  with 
experiment  are  the  overpressure -distance  curve  for  the  main  shock,  and 
the  radius-time  paths  of  the  main  shock,  second  shock,  and  contact  surface. 

Figure  4  presents  a  plot  of  the  calculated  overpressure -distance 
curve  for  TNT  as  a  solid  curve.  Individual  points  are  plotted  for  the 
experimental  data.  (Two  types  of  data  reduction  were  done  by  Fisher  in 
order  to  obtain  the  peak  pressures :  pressures  were  obtained  from  velocity¬ 
line  data,  and  also  from  pressure-time  records.  Both  types  of  measure¬ 
ments  are  recorded  in  Figure  4.) 

The  general  behavior  is  similar  to  that  of  the  previous  pentolite 
curves;  that  is,  beyond  a  certain  point  the  calculated  overpressure- 
distance  curve  remains  fairly  consistently  about  10  per  cent  below  the 
experimental  values. 

Figure  5  is  a  radius -time  plot  for  TNT,  presenting  as  solid  curves 
the  calculated  paths  of  the  main  shock,  the  second  shock  (after  it  has 
emerged  from  the  detonation  products)  and  the  contact  surface.  The 
general  trend  of  the  experimental  data  for  the  main  shock  suggests  a 
behavior  similar  to  that  of  pentolite;  that  is,  at  early  times  the  experi¬ 
mental  points  lie  below  the  theoretical,  crossing  over  eventually  and 
remaining  slightly  above  the  calculated  curve. 

A  new  feature  of  these  curves ,  not  available  for  the  pentolite 
results ,  is  the  capability  of  comparing  with  experiment  the  space-time 
paths  for  the  second  shock  and  for  the  contact  surface.  Figure  5  shows 
that  the  experimental  second  shock  path  lies  about  40  per  cent  above 
the  calculated  one,  while  the  experimental  contact  surface  path  is  about 
a  factor  of  two  greater  in  radius  than  the  calculated  one.  These  are 
serious  deviations,  of  a  much  greater  magnitude  than  the  discrepancies 
in  the  behavior  of  the  main  shock,  and  their  removal  may  require  some 
fundamental  changes  in  the  calculation  and  in  the  model  vised. 

*  Since  Rudlin's  measurements  of  the  contact  surface  are  unpublished,  we 
mention  here  that  they  are  photographic  observations  of  the  luminous 
front,  made  on  high-speed  Ektachrome  color  film. 
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DISCUSSION 

Several  possibilities  can  be  considered  in  an  attempt  to  explain  the 
various  discrepancies  noted  between  the  experimental  and  theoretical 
results.  One  which  is  easily  disposed  of  involves  the  difference  between 
the  loading  densities  used  for  the  experiments  and  for  the  theoretical 
calculation.  Various  experimenters  used  densities  ranging  from  q,  * 

I.5I  -  1.625  gm/cc  for  TNT,  while  the  theoretical  calculation  used  q,  * 
I.625  gm/cc.  A  numerical  calculation  for  q,  “  1-57  gm/cc  showed  essen¬ 
tially  no  difference  with  the  q,  -  I.625  gm/cc  calculation.  Furthermore, 
examination  of  the  experimental  points  on  Figures  4  and  5  fails  to  dis¬ 
close  any  dependence  on  initial  density,  so  that  it  is  probably  safe  to 
discount  this  particular  explanation. 

It  may  also  be  argued  that  a  change  in  the  equation  of  state  of  the 
detonation  products  may  relieve  the  situation  somewhat.  To  test  the 
effect  of  such  a  change,  a  calculation  has  been  carried  out  in  which 
the  LSZK  equation  of  state  has  been  replaced  by  an  ideal  gas  equation 
of  state,  while  the  total  energy  of  the  explosive  (in  this  case,  a  1-lb 
sphere  of  pentolite)  is  kept  constant.  The  results  are  presented  in 
Figure  6,  where  the  solid  curve  represents  the  LSZK  calculation,  and  the 
dashed  curve  represents  the  ideal  gas  calculation. 

The  ideal  gas  calculation  results  in  a  main  shock  which  is  essen¬ 
tially  coincident  with  the  LSZK  main  shock,  and  a  second  shock  whose  path 
lies  about  10  per  cent  above  the  path  of  the  ISZK  second  shock.  Because 
of  the  fact  that  the  second  shock  leaves  the  detonation  products  earlier 
in  the  ideal  gas  calculation,  the  path  of  the  contact  surface  for  later 
times  is  shifted  to  the  left.  However,  the  maximum  radius  attained  by 
the  contact  surface  remains  unchanged,  as  does  its  final  asymptotic 
position.  From  this  single  numerical  experiment  we  can  thus  draw  the 
tentative  conclusion  that  the  main  shock  behavior  and  the  maximum  radius 
attained  by  the  contact  surface  depend  only  on  the  initial  energy  of  the 
explosion.  The  second  shock  position,  however,  does  depend  on  the  equa¬ 
tion  of  state  chosen  for  the  detonation  products.  Nevertheless,  although 
it  does  seem  possible  to  affect  the  second  shock  behavior  by  manipulating 
the  equation  of  state  of  the  detonation  products,  it  would  probably  require 
quite  a  radical  revision  to  achieve  the  corrections  necessary  to  match 
experimental  results .  Furthermore ,  as  long  as  the  total  energy  of  the 
explosion  remains  the  same,  it  is  doubtful  if  the  discrepancies  in  the 
main  shock  and  contact  surface  can  be  removed  by  this  procedure.  Since 
there  exists  some  disagreement  on  precise  values  for  the  heats  of  detona¬ 
tion  of  explosives,  one  might  try  to  correct  the  main  shock  behavior  by 
increasing  the  heat  of  detonation  by  a  suitable  amount  in  the  theoretical 
calculation.  This  does  result  in  a  suitably  adjusted  main  shock  (at 
least  in  the  region  where  the  main  shock  pressures  were  initially  too 
low)  but  shifts  the  position  of  the  contact  surface  and  second  shock  by 
approximately  the  same  percentage  as  the  main  shock,  which  is  hardly 
enough  to  account  for  the  discrepancies. 
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Perhaps  a  reasonable  approach  is  to  combine  the  last  two  ideas  by 
adjusting  the  heat  of  detonation  to  give  the  correct  main  shock  behavior, 
and  to  try  to  find  a  suitable  equation  of  state  for  the  detonation  products 
which  will  remove  the  remaining  discrepancies  in  the  second  shock  and 
contact  surface.  On  the  other  hand,  tie  realization  that  the  second  shock 
behavior  is  affected  strongly  by  the  state  of  the  detonation  products 
leads  one  to  remember  that  the  theoretical  calculation  is  a  highly  idealized 
one  which  omits  consideration  of  many  phenomena.  In  particular,  the 
numerical  calculation  does  not  take  into  account  possible  reaction-zone 
effects,  turbulence,  jetting,  incomplete  combustion,  mixing,  afterburning, 
etc.  It  is  quite  conceivable  that  some  or  all  of  these  effects  are  impor¬ 
tant  and  that  a  more  accurate  prediction  can  be  obtained  only  when  a  more 
sophisticated  model  is  postulated  for  the  behavior  of  the  expanding  gases. 
More  detailed  experimental  determination  of  the  explosion  phenomena  may  be 
required  for  the  creation  of  a  satisfactory  theoretical  model. 
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FIG.  1  PEAK  SHOCK  OVERPRESSURE  VS.  RADIAL  DISTANCE,  1-LB.  SPHERE  OF  PENTOLITE 

AT  SEA  LEVEL  CONDITIONS 
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FIG.  2  RADIUS-TIME  CURVES  FOR  1-LB  SPHERE  OF  PENTOLITE  AT  SEA  LEVEL  CONDITIONS 
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FIG.  3  POSITIVE  DURATION  VS.  RADIAL  DISTANCE,  1-LB.  SPHERE  OF  PENTOLITE 

AT  SEA  LEVEL  CONDITIONS 
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FIG.  4  PEAK  SHOCK  OVERPRESSURE  VS.  RADIAL  DISTANCE,  1-LB.  SPHERE  OF  TNT 

AT  SEA  LEVEL  CONDITIONS 
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FIG.  5  RADIUS-TIME  CURVES  FOR  1-LB.  SPHERE  OF  TNT  AT  SEA  LEVEL  CONDITIONS 
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FIG.  I  PEAK  SHOCK  OVERPRESSURE  VS.  RADIAL  DISTANCE,  1-LB.  SPHERE  OF  PENTOLITE 

AT  SEA  LEVEL  CONDITIONS 
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FIG.  3  POSITIVE  DURATION  VS.  RADIAL  DISTANCE,  1-LB.  SPHERE  OF  PENTOLITE 
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FIG.  4  PEAK  SHOCK  OVERPRESSURE  VS.  RADIAL  DISTANCE,  l-LB.  SPHERE  OF  TNT 

AT  SEA  LEVEL  CONDITIONS 
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